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ABSTRACT 


The Autonomous Underwater Vehicle (AUV) must be able to operate in various 
shallow water sea state conditions. In order to have a precise navigation and steering 
system, and efficiently place charges on underwater mines, the AUV must be able to sense 
and overcome hydrodynamic forces which are caused by waves. This thesis establishes a 
model of sea state conditions based on spectral analysis, and uses the model to predict 
future knowledge of the sea. This prediction is determined by the random white noise 
output of a discrete, digital filter. The development of the discrete, digital filter is 
described herein. The Pierson-Moskowitz (P-M) spectrum which models seaway 
elevations using linear wave theory is used as a target spectrum which the filter will track. 
Cross-correlation between the P-M target spectrum and digital filter have shown that a 
reasonably accurate estimate of wave elevations can be predicted one full wave period into 


the future. 
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I. INTRODUCTION 


A. HISTORY ON MINE APPROXIMATION COUNTERMEASURES 


Once the Navy’s primary mission shifted from “Sea Control” (blue water 
operations) to “From the Sea” (shallow water operations), the increased requirement for 
amphibious warfare operations demanded that increased interest be placed in the area of 
mine countermeasures - specifically, mine warfare in shallow water. Because of the 
immense level of tension and somewhat frequent conflicts recently arising in the Middle 
East, the interest of the United States security has often been channeled toward the mine 
infested waters of the Arabian Gulf. Over the past ten years or so, the United States has 
discovered and cleared a large number of underwater mines in this area. Furthermore, 
there have been some adverse developments occur as an aftermath of ships encountering 
explosive mines. To this end, we must continuously remain attentive of this phenomenal 
threat to our forces and continue , to seek improvement toward measures of combating 
this hostile and unwavering enemy by devising safe and efficient tools to prevent mines 
from presenting future problems, and hindering the effective operation of Amphibious 
Warfare. 

As a measure to alleviate the danger of underwater mines, the Naval Postgraduate 
School Center for Autonomous Underwater Vehicle (AUV) Research through 
Mechanical, Electrical Engineering, and Computer Science Departments are continuously 


researching new methods to overcome this ongoing challenge. 
B. MINEHUNTING IN SHALLOW WATER (10 TO 40 FEET) 


With the new mission (From the Sea) in place, today’s Modern Navy, more than 
ever before has focused on means to intelligently and efficiently counter the challenge 
presented by floating mines by employing countermeasures in shallow water. Mine 


hunting in shallow water has presented continuous challenges to the safe and effective 








operation of the United States Naval Forces. For maximum flexibility the mine 
countermeasure efforts should be covert, cost effective and relatively quick, (Richwine, 
1995). The countermeasure methods currently used by Naval Forces include the 
involvement of humankind, noisy platforms such as ships and helicopters. These 
aforementioned methods require a high risk of human sacrifice, and compromises 
covertness respectively. As another measure, marine mammal systems and special military 
forces are capable of operating covertly. However, they are scarce resources that require 
extensive training pipelines, (Hunt, 1995). 

Owing to the robust advancement in technology, countermeasures using 
autonomous systems are currently being built and tested at the Naval Postgraduate 
School. In previous studies at the Navy Postgraduate school, an autonomous underwater 
vehicle (Phoenix) has been developed and operated in a test tanks where sea state 
conditions were negligible. In order to advance the studies and operability of the AUV, 


wave and current forces must be considered and incorporated in the test environment. 


C. DISTURBANCE REJECTION 


The main source of the dynamic forces encountered by underwater vehicles or 
submerged vessels are wave and current induced. These disturbance forces arise from 
buoyant and inertial effects due to ocean wave kinematics. In the past, most sea state 
models were predicted using an erogodic process. To best assess the performance of an 
AUV in shallow water, it is necessary to generate a dynamic model of the velocity and 
acceleration particles of water to model the wave forces. Although frequency 
decomposition is a well known method to simulate time history variation of the ocean 
surface due to waves, the digital representation of the sea surface elevation and pressure 
profile process was thought to offer advantages, and is referred to as a time series of wave 
history. Because there is a strong uncertainty of random ocean wave fields, in the past, a 


spectral approach has been best suited for describing the intensity and frequency 


characteristics of ocean wave fields. 














Recent research on simulating ocean disturbance was conducted by Spanos during 
the an Offshore Structures Conference, (Spanos and Mignolet, 1986). His research 
expanded on the use of the Linear Prediction Theory using Autoregressive (AR) and 
Autoregressive Moving Average (ARMA) Models. With the former, standard method, 
considerable numerical difficulties were encountered in getting realistic simulations, while 
the latter methods experienced difficulties accurately approximating the target (P-M) 
spectrum, and an accurate approximation could only be feasible when using a high order 
AR scheme, (ibid). Spanos used dimensionless coefficients in his version of the ARMA 
model in order to utilize the filter for multiple wind conditions and sea states. 

It was found that the spectrum of the time series generated by the AR model 
exhibited considerable narrowband fluctuations with respect to the P-M target spectrum. 
The fluctuations are thought to be associated with pole locations of the transfer function 
with respect to the unit circle. The dimensionless coefficients used by Spanos resulted in 
high frequency models that fail to match the P-M spectrum. Additionally, the form of 
dimensionless coefficients did not properly illustrate the behavior of ocean waves under 
different sea state conditions, compared to the nominal one chosen for the normalization. 
As it turns out, the mathematical form of the P-M spectrum cannot be realized exactly in a 


linear filter transfer function. 
D. OBJECTIVE OF THESIS 


Even though the spectral analysis of ocean wave-induced motions have been 
studied over a period of recent years, problems and limitations continue to exist. These 
limitations have made it difficult to predict sea state conditions; and consequently have 
made the task of efficiently operating Autonomous Underwater Vehicles challenging. 

For an Autonomous Underwater Vehicle (AUV) to operate in shallow water with 
adequate reliability, wave forces and their effects on the AUV should be modeled 
mathematically with an adequate degree of accuracy. A major component of the 


mathematical model of these loads/forces experienced by underwater vehicles involves 








ocean wave kinematics, (Sarpkaya and Isaacson, 1981). In fact, over a period of years 
the Pierson-Moskowitz (P-M) spectrum has dominated the statistical description of ocean 
wave kinematics, (Pierson and Moskowitz, 1964). Thus, in analyzing the impact of these 
forces on underwater vehicles, it is worthwhile to generate a time record of values of wave 
kinematics that is compatible with the P-M spectrum. Traditionallythis has been achieved 
by relying on the harmonic superposition method, (Hudspeth and Borgman, 1979, 
Shinozuka, 1972, Shinozuka and Wai, 1979), which can be cumbersome and costly in 
terms of computation time. Recently, with emphasis placed on improving the 
computational time, Spanos developed a procedure based on digital filtering techniques. 
His method involved the design of an autoregressive filter (AR) and autoregressive 
moving average (ARMA) filter that consisted a transfer function whose squared modulus 
was a close of the target P-M spectrum, (Spanos and Hansen 1981, Spanos 1983). His 
digital filtering technique demonstrated promise in computational efficiency, nevertheless, 
it proved to be a delicate task for spectral matching, (Spanos and Mignolet, 1986). 

In this thesis, seaway dynamics will be modeled as a linear system with a random 
white noise input, where the linear system transfer function 1s approximately mapped to 
the P-M spectrum. Further, this work will use a digital filtering technique from two 
transfer functions to approximate the P-M spectrum. Initially, the effectiveness of 
Spanos’ AR and ARMA methods will be determined by testing the stability through 
eigenvalue analysis. These eigenvalues are to be modified as necessary to achieve required 
stability. This model will then be compared to the output of a digital filter generated from 
the transfer function of an eighth-order system. Once the results of the eight-order 
system are compared to the results of the Spanos model for best mapping of the P-M and 
pressure profile spectra, at specified water depths, the system (eighth-order or Spanos) 
with the most accurate approximation of the P-M target spectrum will be further studied | 
in a filter system and used to predict sea states and pressure responses at future times 


based on current and previous measurements. The accuracy will be measured using cross- 


correlation between the target spectrum and the dynamic equation/innovator. 








II. LINEAR PREDICTION THEORY (LPT) 


A. HISTORY AND FORMULATION OF LPT 


Digital representation of sea surface elevation and pressure profile processes is 
called a time series. The probabilistic description of ocean and wave induced responses 
was introduced during the 1950s, (Spanos and Mignolet, 1986). Prior to the introduction 
of the Linear Prediction Theory (LPT), the traditional techniques for generation of 
simulated digital sea wave height, velocity, and acceleration records were based primarily 
upon the principle of wave superposition. For each point in time domain, this technique 
generally required the computation of a sum of harmonic waves, (Spanos and Hansen, 
1981). Nevertheless, because of successful research efforts in the areas of simulation and 
processing of digital data in the time domain, modeling of sea wave time history with the 
Linear Prediction Theory has proven to be quite accurate. The LPT method provides an 
efficient procedure of generating a collection of records possessing a power spectrum 
which approximates a known target spectrum, (ibid). In the LPT, the records of sea 
waves are approximately compatible with the Pierson-Moskowitz (P-M) target spectrum, 
and are obtained as the output of a recursive digital filter with a random white noise 
input. The problem of base period repetition of wave records is thus avoided. 
Essentially, this technique simply realizes the approximation of the Pierson-Moskowitz 
(P-M) spectrum as the output of a digital filter to white noise excitation. More details on 
the development of the P-M Spectrum formulation will be addressed later. 

When modeling the spectrum, the velocities and acceleration particles for sea 
wave simulations can be simulated utilizing the Linear Prediction Theory and the 
following relationships for the horizontal velocity P,(@) and acceleration P,(q@) spectra, 


(Shinozuka, 1972) 


P,(w) = @|G(a,z, H)| P(o) 


v 











and 


P,(o) = w'[G(a, z, H)| P(o) 


where Pq) is the surface elevation spectrum, and G(a,z,H) is a function of @ depending 
as well upon the total water depth (/7) and the vertical position (z) above the sea floor 
where the velocity and acceleration simulations are desired. The independent variable @ 
represents frequency in radians per second. Using Py(w) and P,(@) as target spectra, LPT 
can be applied to obtain time domain velocity and acceleration records. This method of 
simulation is based on an all-pole recursive filter, the input, as stated above, of which is a 
random white noise process. A one-time computation of the digital filter difference 
equation coefficients are required. Once these coefficients are computed for a particular 
target spectrum, they can be stored for future usage as necessary. The generation of the 
output time domain record requires only simulation of white noise digital samples and 
application of the difference equation. No computations of any trigonometric or other 
functions are required after the initial coefficient generation. 

The LPT was utilized to digitally simulate the time domain sea wave height 
records. The target P-M spectrum was treated as one-sided extending over the frequency 
range (0,@;). White noise records were numerically simulated and used as input to the 
difference equation representing the digital filter. The difference equations and methods 
for determining the difference equation coefficients can be cited in (Spanos and Hansen, 
1981). The Fast Fourier Transform (FFT) of the digital time record generated by the 
difference equation was used to calculate the response spectrum of the digital filter 
output. 

The following figures represent the results obtained using LPT. Figure | 
illustrates an example record for an upper frequency @» = n/2 rad/sec, filter of order p = 
25, and wind velocity of 40 knots. Figure 2 illustrates a portion of a record which is the 
output of a digital filter of order 50, upper frequency of @ = 7/2 rad/sec, and wind 


velocity of 40 knots. Figure 3 is the average power spectrum of 300 records of 1000 time 








increments length for an upper frequency @, = 1/2 rad/sec and wind velocity of 40 knots. 
A sampling interval of T = 2 seconds between pulses of the filter output to white noise 
input was used. This procedure was repeated for simulations of lower order filters as 
well. Noticeably, a trend of decreased mapping between the P-M target spectrum and 
power spectrum waves of the filter output was obvious. Figures 4 through 7 illustrate 
records and power spectra of outputs of filters of orders 25 and 40 respectively, with 
frequency @,= 1/3 rad/sec, a sampling interval of T = 3 seconds, and a wind velocity of 


40 knots. 
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Figure 1. Normalized autocorrelation, (Spanos and Hansen, 1981). 
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Figure 2. Portion of a LPT-simulated sea wave height record, 
(Spanos and Hansen, 1981). 
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Figure 4. Portions of a LPT-simulated sea wave height record, 
(Spanos and Hansen, 1981). 
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Figure 5. Power Spectrum of LPT-simulated sea wave height record, 
(Spanos and Hansen, 1981). 
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Figure 6. Portions of a LPT-simulated sea wave height record, 
(Spanos and Hansen, 1981). 
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B. MODELING OF THE SYSTEM AND REAL TIME-FORWARD 
PREDICTION 
The general block diagram model of the proposed overall system including the 
target spectrum transfer and the transfer function of the dynamic equation/filter with the 


innovator gains illustrated is depicted below in Figure 8. 





Figure 8. Block diagram of system 


where w/(s) is the white noise input, Y (s) is the filtered estimate of the under wave 


pressure response, Y(s) is the measured pressure, E,(s) is the error between Y(s) and y (s) ; 


and K’s are the system gains. The independent variable s represents the Laplace variable. 


If the following expression is stable then the error (¢) between the filter estimate and the 


wave response is a bounded value as foo, 


G(s) 


1+ KG(s) 











Y(s), wave motion (i.e. measured pressure signal) is assumed to arise from a random white 


noise driven system with a dynamic transfer function where, 


|G) = S@) 


Spanos developed a transfer function (in terms of z-transformation) which is a simple 
version of the more general higher order general Autoregressive Moving Average 


(ARMA) model. 


> b a 


G(z)= = - 
1+ > a,z* 
k=) 


This transfer function in terms of z-transformation notation was used as the output of the 
digital filter which is represented by the symbol Y(0). 

In Spanos’ dimensionless model the coefficients a; and b; were optimally chosen to 
minimize error between the Pierson Moskowitz target spectrum equation and the transfer 
function of the ARMA filter. Spanos presented the coefficients in dimensionless forms in 
order to use them in various structural dynamic applications. Nevertheless, when 
attempting to model a dynamic equation for the target P-M spectrum, the transfer function 
generated by the coefficients of the dimensionless model resulted in a narrow band, high 
frequency spectrum that failed to adequately resemble and map the P-M spectrum. The 
high frequency results are illustrated in Figure 9. Additional problems were encountered 
with normalization - it was difficult to reconstitute and stabilize the system, and the 


transfer function had unstable zeros, more details on Spanos’ ARMA in following 


chapters. 
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Mathematical calculations determining the error between the P-M target spectrum 
(wave and pressure responses) and the dynamic equation/filter have been performed to 
verify the validity of the system’s block diagram. Figure 10 illustrates a modified block 
diagram of the system which emphasizes how the error (é) is measured schematically. 
Assuming that the seaway can be represented as the output from a white noise driven 
filter, generally, the error (€) goes to zero when the wave and pressure responses G(s) 


equals or matches the response of the dynamic equation S(q@), or simply : 


s—>0; |c(o)| = S(a) 
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Figure 9. Time Record of Spanos’ Dimensionless ARMA Model 
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- Figure 10. Block diagram of Forward Prediction System, where s = seaway internal 
dynamics and P,,, = ||5| 
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Assuming, for the moment that all white noise inputs are unity, 
e=Y-Y=-S*w+G(w-Ke) 


therefore, 


(6-3) 
(1+ GK) 


e(1+GK)=(G-S)w > &é&= 
From examining the numerator of the above expression, it is obvious that when Gai the 


numerator is zero, therefore e=0, and if GK is stable, ¢ is bounded and depends on 
IG - S| <|[G]-|s| 
since G and S are both positive, real. 


The equation to follow is given in terms of a transfer function, and has been used to 
mathematically model the transfer function of the dynamic system. In the work of Spanos, 
the equation is given as the digital filter of the following transfer function represented in 


terms of z-transformation: 


> bz 
Gq) = —=,— 
1+ De a,z-* 


k=] 


In this domain, z is the complex variable z = e”™" and s is the continuous time Laplace 
| P p 


variable. For frequency inputs, s = j@ (q is the radian frequency), in terms of the 
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corresponding frequency response of the discrete time system, we substitute z = é/" and 


the magnitude of the system transfer function becomes 




















which depends on both, @ and T. 


For Time Domain Representation: 


Take G from the model and convert it to the following equation in the continuous time 


domain. 

X(t) = Ax(t)+ Bu(t), y = Cx 
where 0 < t < «, and x(t) is a continuous state and w(f) is a continuous input of the 
interval 0 < t < co, with the initial values x(t = 0), u(t = 0) and u(t > 0) being a series of 


zero order hold data such that at t, = kdt, u(t,) is a data sequence: Between the intervals 


ty <t< tys1, where u(t) =u, (Figure 11). 
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Figure 11. Represents Discrete Time Domain 


The solution for x(t=t,) is given by the following integral 


{ x(1)dt = { Ax(t)dt + { Bul at 
ty t 


ty 
giving 
Xuer= (C7 )xn + (2 - Ya" bu, => Xe = Gt 
where 


g=e" and y= (e" -1 )A"B. 


This converts x = Ax + Bu to the following system in discrete sampled time domain 
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x(t) = Ox(f -1)+ Tut) 
y(t -1) = Cx(t -1) 


Initially, Spanos used pole placement to design the control matrix, K . However, 
the Linear Quadratic Regulator (LQR) proved to be more convenient for determining the 
values of the gain matrix, leading to the estimator ( u in this case is zero mean, white 


noise) 


Xy4 = PX, +K(p, ~ p); p= Cx, 


where p; is defined as the measured pressure from the Pierson-Moskowitz and pressure 


profile spectrum. 
Simplifying the above equation yields: 
£,., =(®- KC)z, + Kp,; p, = CX, 
The best estimate/prediction forward in time is then using the above model without 


data correction obtained through the following equation from the discrete time domain. 


This forward prediction is always based on the estimated quantity of £(t) : 


P(t +k) = 2*Cx(t) 
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C. PROBLEMS ENCOUNTERED WHEN USING LINEAR PREDICTION 
THEORY (LPT) 


This method superseded the harmonic superposition method. However, it proved 
to be quite cumbersome for data processing machines requiring immense computational 
time. Problems were noted with the fluctuation of the LPT-generated spectrum from the 
target spectrum. Figures 3,5 and 7 illustrate these fluctuations. The LPT is not always 
computational efficient for simultaneous simulation of sea wave velocity and acceleration 
records at several times, and points with any combination of horizontal and vertical 
coordinates, (Spanos and Hansen, 1981). Considering the behavior of the power spectra 
when mapping it with the target spectrum as the order of the filter p was decreased, it 
was discovered that increasing the filter order p, the local frequency breadth of the 
spectral fluctuation decreases, (ibid). 

Using an approximation of the P-M spectrum as the output of digital filters to 
white noise excitation has been proven to have numerical difficulties when approximating 
it by the Autoregressive (AR) models, (Spanos, 1983). However, when using the AR 
model as an initial spectral approximation, it proves to yield an efficient autoregressive 


moving average (ARMA) model, (Spanos and Mignolet, 1986). 
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mW. APPLICATION TO WAVE MOTION PREDICTION 


A. PIERSON-MOSKOWITZ (P-M) SPECTRUM 


A semi-empirical expression for the frequency spectrum of fully developed seas is 


given by, 


so) oa ie) 


This expression was proposed by Pierson and Moskowitz for the height of fully 
developed seas, (Pierson and Moskowitz, 1964). The a@ and # coefficients are 
nondimensional parameters defining the spectrum. The P-M spectrum was developed for 
wind velocity at a standard height of 19.5 meters above the free surface. This two 
parameter spectrum is sufficiently general to fit quite a few observations and is consistent 
with theoretical predictions of the high frequency limit. Where most common values for 
the nondimensional parameters are 

a. = 8.1x 10° 
B=-0.74 
The gravity constant is 
g = 32.174 ft/sec” 
The two variables in the expression are 
@ = frequency (rad/sec) 


V = wind speed (ft/sec) 


The Pierson and Moskowitz expression can be further reduced by considering the 
relationship between the significant wave height, Hz, and wind speed, V. The following 
expression represents the aforementioned relationship: 
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~ 1/2 2 
H,, = ad s(0)do) = 02092— 
0 


Substituting this relation, along with the specified nondimensional parameters and gravity 
constant, into the Pierson and Moskowitz spectra equation yields a simplified version of 


the P-M expression which was used to generate the P-M and pressure profiles for the 


target spectrum (S(@)). The simplified expression is 


S(@) 5 HE? xo’ 


81x107 x 32174" | pa x al 
= —_______—_——_- x exp| —0.74 ———_.—_,-_ 
@ 


where H was for (the mean water depth) shallow water, linear waves for this particular 
model is 7 meters (approximately 23 feet). (Papoulias, 1993) 

Typical plots of varying significant wave heights (Hs) versus frequency and period 
for P-M spectra are illustrated in Figures 12 and 13 respectively. 


B. _PIERSON-MOSKOWITZ (P-M) AND PRESSURE PROFILE TARGET 
SPECTRUM 


In performing digital time series analysis, the sampling period 7 and the cut-off 
frequency @, are related through the Nyquist relation T= wa. In order to generate the 
appropriate target spectrum at a specified operating depth, in this case 5.2 meters 
(approximately 17 feet) below the surface, the pressure at this depth had to be considered. 
To this end, the P-M expression was multiplied by G(@,z,H) which resulted from the linear 
wave theory. Plots consisting of varying significant wave heights (Hs) versus frequency 
and period for P-M and pressure profile spectra are illustrated in Figures 14 and 15 


respectively. 
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Figure 12. Pierson-Moskowitz (P-M) spectra for significant height (H,) of 8, 6, 4, and 3 
feet versus frequency. 
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Figure 13. | Pierson-Moskowitz (P-M) spectra for significant height (H,) of 8, 6, 4, and 3 
feet versus periods. 


28 








. 
. 


—-- Hs =8 feet 

* Hs ='6 feet 
| } 7 : / |- — Hs =4 feet 
3.5 3 e 3 .. | ——.Hs. ='3 feet 


OO 


NO 


—k 
O1 


S(w) mean square feet per rad/sec 
NM 
O1 





Frequency—w (rad/sec) 


Figure 14. | Pierson-Moskowitz (P-M) and pressure profile spectra for significant height 
(H,) of 8, 6, 4, and 3 feet versus frequency. 
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Pierson-Moskowitz (P-M) and pressure profile spectra for significant height 
(H,) of 8, 6, 4, and 3 feet versus periods. 
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22 


The Depth Modifier for pressure variation expression is 


cosh(211(H — eA 
Pwave(t,X,Z)=PZZtNapg ee a sin an( 4 - *) 
cosh( 211 ) roa 


Revised Depth Modifier for the pressure variation expression neglecting gravity 


cosh(2II(H ~ 2)/ 
G(z,H,A) = ~ cosh(HHZ) 


The expression for generating P-M and pressure profile target spectrum is 


cosh(2I1(H — zy 


‘ cosh 21H, | 


: H? xo* 


81x10? x 32174 (22174 x oan) 
ne x exp] -0.74| 
@ 





Taking the regular wave theory, and assuming a horizontal sea bottom and a 


free-surface of infinite horizontal extent, a linear wave theory for propagating waves can 


be derived. This derivation applies the free-surface and sea bottom conditions, (Faltinsen, 


1990), and the Laplace equation V’*¢ = 0. The mean water depth will be defined by the 


constant H. The following expression then gives the connection between wavelength A 


and wave period 7 for a finite water depth. 


pn See rte eal 2 
2 A 


3] 





Substituting @ = 27/A in the above expression and simplifying the nonlinear quadratic 


equation yields 


f= (28 - gtanh(201H7/ ) =0 


where all values of lambda A were found as functions of frequency @, for this particular 
depth condition see Table 1. 

For this particular research, the wavelength and frequency connection equation 
was solved for a specified mean water depth (H) of 7 meters. The quadratic was solved 
using the recurrence equation solver in conjunction with the floating-point arithmetic 
solver from the Maple Software Application. The frequency range for lambda Aas a 
function of @ is 0.2 < @(i) <3, in increments of 0.05. 

Thorough details regarding the derivation of linear wave theory for propagating 
waves are beyond the scope of this work, however, information can be found in many 
textbooks in fluid mechanics (for instance Newman, 1977, chapter 6), (Faltinsen, 1990). 
The equations for free-surface and sea bottom conditions, along with the relations of 


wavelength / and period 7 can be cited in (Faltinsen, 1990). 
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Table 1. Lambda (4) defined as function of omega (@) for water depth of H = 7 meters 
(approximately 23 feet). 
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C. SPANOS’ DIMENSIONLESS ARMA FILTER 


Recently, Spanos developed a digital filtering techniques that involved the design 
of an autoregressive moving average (ARMA) filter having a transfer function whose 
squared modulus was assumed to be a close approximation of the target P-M spectrum. 
The transfer function of the ARMA(/,m) process y is a discrete random process whose n™ 


sample can be obtained from the previous one using the following expression: 


I m 
Yn = —DL ans - yb nj 
=] j=0 
where @ denotes a bandlimited /-a@,@] white noise process with autocorrelation 


function, (Spanos and Mignolet, 1986). 


Spanos also used the process y as the output of a digital filter of a transfer function in the 


terms of z-transform notation. The input to this transfer is a discrete band limited white 


noise process @. 


bz? 
H (z ) = = 
1+ > a,z* 

k=1 
This transfer function was examined using the dimensionless coefficients in Table 
2. After comparing the power spectrum of the transfer function and developing a wave 
record, it was discovered that this process resulted in a very high frequency record which 
was far from mimicking the behavior of ocean waves. This behavior was cited earlier in 
Figure 9. Consequently, another attempt to utilize this transfer function for matching the 


P-M Spectrum was made by modifying the numerator in order to best represent wave 


motions. 
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Manipulating the latter transfer function by taking the z transforms to get a z-domain 


(ARMA) model of system dynamics yields 
zX=oX+yYU => (z-o)X=yU 


where ¢=e%7 €R™ and y © R™ 


therefore, 
-1 


o@)=F=[-o'7 => OG)=7 Es 


where y and ¢ are m and n"-order polynomial expressions. 


Spanos ARMA model transfer function is simplified to a more general higher order 


ARMA model expression such as: 


x bz? 
G(z)=—= 


— a 
1+¥a2" 
k=) 


in which b; and a, are related to yand ¢@. 
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Table 2. Spanos’ P-M Wave Elevation Spectrum ARMA Coefficients (dimensionless). 


-1.2058x10" 2.1836x10" 
-2 
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D. EIGHTH-ORDER FILTER 


Taking into account Spanos used an eighth-order transfer function to generate the 
system of dimensionless coefficients, a eight-order filter was assumed to be best suited for 
generating a filter to model the target spectrum in this work. However, spectral densities 
or wave spectra for transfer functions of different orders, ranging from second to eight- 
order systems were evaluated and compared to the target wave spectrum. After 
evaluating the spectrum of all aforementioned ordered systems, the spectrum of the 
eighth-order system, as expected, most closely reflected that of the target wave specrtum. 
The comparison between the target wave spectrum and the eighth-order system’s 


spectrum is depicted in Figure 16. 


The following expression was used to formulate the eighth-order transfer function for the 


dynamic equation which was used as filter to forward predict seaway surface elevation. 


+20 % )+(%,) 


Generally, the eighth-order filter process P is a discrete random process whose Ke” 


sample can be obtained from the previous sample k in the following way 
Xp = Py +q;P, = Cx, 


where the state matrix ¢ is determined by the transfer function and q is a random white 


noise input with zero mean. 
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Figure 16. Predictor Spectrum versus P-M Spectrum. 
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The filter innovator (gain matrix K) for this system was generated by using the 
LOR for determining the stability of the system. The frequency was chosen according to 
the peak frequency of the target wave spectrum. This value was determined to be 
approximately 1.31 rad/sec. It is worthwhile to note that the peak frequency for the P-M 
and P-M/pressure profile spectra was approximately 1.2 rad/sec. This frequency value, 
even though the difference is relatively small when comparing to that of the P-M spectrum 
alone, didn’t match the target spectrum as well as the frequency of the P-M spectrum. 

Tests conducted with the dynamic system using the two different peak frequencies 
were conducted. The frequency from the P-M spectrum predicted sea condition out to a 
further time with better accuracy. The filter gains for the LQR system were generated to 
match the target wave spectrum. Noticeable was the small difference between the filter of 
the pressure and target spectrum. 

For this work, significant wave height of 3 feet was used for both the P-M target 
spectrum and pressure profiles. Figures 17 through 20 illustrate the frequency and period 
plots, depicting the significant wave height at this particular depth, for each spectrum. 
Figures 21 and 22 illustrate the relationship between the P-M spectrum pressure profiles 
by normailizing and superimposing them on the same plot. Figures 23 through 25 
represent time records of the P-M spectrum, pressure profiles, and eight-order filter 
spectrum respectively. When comparing figure 25 to figure 9 (Spanos dimensionless time 
record), it is obvious that they differ immensely. The a; and b; coefficients of the eighth- 


order filter are listed in Table 3. 
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Table 3. Eighth-Order Filter Coefficients. 


40 








0.7 


2 
op) 


Oo 
on 


© 
SN 


© 
oF) 


0.2 


S(w) mean square feet per rad/sec 


0.1 





0 0.5 1 1.5 2 2.5 3 
Frequency—w (rad/sec) 


Figure 17. Frequency plot of P-M Wave Spectrum for significant wave height (H 1) 
of 3 feet. 
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Figure 18. Period plot of P-M Wave Spectrum for significant wave height (H 1) of 3 
feet. 
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Figure 19. Frequency plot of Pressure Profile Spectrum for significant wave height 
(H;;3) of 3 feet. 
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Figure 20. Period plot of Pressure Profile Spectrum for significant wave height (Hy) 
of 3 feet. 
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Frequency plot P-M spectrum superimposed on Pressure Profile Spectrum 
for significant wave height (Hi;) of 3 feet (normalized spectra). 
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Figure 22. Period plot P-M spectrum superimposed on Pressure Profile Spectrum for 
significant wave height (Hi) of 3 feet (normalized spectra). 
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Time History of Pierson-Moskowitz Spectrum reflecting random output 
for significant wave height (Hi, ) of 3 feet. 
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Figure 24. Time History of Pressure Profiles reflecting random output for significant 
wave height (Hj); ) of 3 feet. 
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Figure 25. Time History of Measured Pressue Profile reflecting random output of 
eighth-order, digital filter for significant wave height (His ) of 3 feet. 
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IV. RESULTS 
To aid in the discussion of results, it is helpful to give the following, 
Definition: 
p(t,(t+7)): estimate of forward pressure at time, (f+ 7) predicted at current time, (1). 


y(t): measured pressure at time, (1). 
y(t +): measured pressure to be found from seastate at time,(t + T). 


Prediction error, (,) at time, (t + t) based on filter data up to time (t): 


e(t, (¢+ t)) = Alt, (t+ t)) — y(t+T) 


Figure 26 illustrates how this error is defined graphically while Figures 27 through 38 
graphically reflect time history plots of the results obtained from the above definitions. 

A measure of effectiveness is the forward prediction using the correlation 
coefficient between digital output of the P-M/pressure profile target spectrum and transfer 
function of eighth-order system. This measure of effectiveness predicted wave and 
pressure responses out to designated number of time steps (number of time steps for each 
measure is reflected on plots of the results). However, in this work, results were taken at 
time, t=5 and t=30 seconds. A fluctuation was observed with the cross-correlation 
results. This fluctuation is mostly believed to be the aftermath of the random phenomenon 
of the seaway time history record which had to be modeled with a random white noise 


output process of a digital filter. As a measure to compensate for the fluctuation and 
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obtain the best results, the output data were calculated fifty iterations and averaged to get 
the final data value. The graph in figure 39 depicts the mean cross-correlation coefficient 
flanked by the standard deviation. 


p(t+7) and p(t+7) were compared over the entire time record, while the following 


expression was used to formulate the cross-correlation coefficient. 


3 p(t +7) p(t,t+ t)/N 


A=" std(p)std(p) 
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Figure 26. | Forward Prediction Diagram 
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Figure 27. Time Record of P-M Target Spectrum (Y) versus Measured Presure (p/, 
(Predictor). 
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Zoomed Plot of Time Record of P-M Target Spectrum (Y) versus 


Figure 28. 
Measured Pressure (p), (Predictor). 
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Time Record of Measured Pressure (p) versus Measured Pressure p(t+ 7). 


Figure 29. 
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Zoomed Plot of Time Record of Measured Pressure (p) versus Measured 


Figure 30. 
Pressure p(t+ 7). 
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Time Record of Measured Pressure p(f+ 7) versus Estimate of Forward 


Figure 31. 
Pressure Bit, (t+ t)) 
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Zoomed Plot of Time Record of Measured Pressure p(f+ 7) versus 
Estimate of Forward Pressure Alt, (¢+ t)) . 
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Figure 33. Time Record of P-M Target Spectrum (Y) versus Measured Pressure (p/, 
(Predictor). 
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Zoomed Plot of Time Record of P-M Target Spectrum (Y) versus 


Figure 34. 
Measured Pressure (p), (Predictor). 
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Figure 35. Time Record of Measured Pressure (p) versus Measured Pressure p(t+ 7). 
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Zoomed Plot of Time Record of Measured Pressure (p) versus Measured 


Figure 36. 
Pressure p(t+ 7). 
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Figure 37. 
Pressure Blt, (t+ r)) | 
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Zoomed Plot of Time Record of Measured Pressure p(t+ 7) versus 
Estimate of Forward Pressure Bit, (¢+ t)) 
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Figure 39. Mean Cross-Correlation Coefficient and Standard Deviation between P-M 
Target Spectrum and Digital Filter. 
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V. CONCLUSION AND RECOMMENDATIONS 


A. CONCLUSION 


In previous studies of hydrodynamic effects on underwater structures, attempts to 
simulate waves using the Pierson Moskowitz spectrum have long been of interest. This 
work has revealed that using an eighth-order filter is a convenient means of matching the 
P-M spectrum and predicting future seaway surface elevation behavior. 

It was discovered that the application of normalized coefficients to model seaway 
for various states proved to be inconsistent with the actual behavior of sea waves. This 
discovery was explored when a comparison of spectral densities and time history plots 
between dimensionless and an eighth-order transfer function coefficients was made. From 
comparing the spectral densities and time history of Spanos’ model and the eighth-order 
filter with the spectral density of the P-M spectrum, it became apparent that Spanos model 
failed to match the P-M spectrum while the eighth-order system proved to be promising , 
matching the P-M spectrum very closely. 

It has been concluded that an eighth-order transfer function provides good 
coefficients for tracking the P-M target spectrum and predicting seaway elevation 
response. What cannot be concluded from this particular work is the degree of optimality 
to which the eighth-order system can actually match the target spectrum. Thus far, based 
on history of previous studies in this area, the work conducted herein has indicated a 
promising direction for better predicting future seaway response. From this research, it is 
indicated that the eighth-order system allows adaptivity and a degree of future prediction 


capability based on past and current data. 
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B. RECOMMENDATIONS 


This filtering system is merely a baseline version for modeling seaway, the gain 
matrix for the innovator was not optimized. That is, the Linear Quadratic Regulator was 
not optimally chosen - the weighting matrices values were based on intuition drawn form 
past experiences without access to experimental data. If the Q and R matrices are chosen 
the predictor would most likely predict further into the future with more accuracy. Also, 
the ability of the prediction model to be adapted to seaway spectra that vary in peak 
frequency and frequency spread should be explored. 

Also, it should noted that the P-M target spectrum is recognized as an asymptotic 
form and it has some limitations, for shallow water applications. Therefore, the 
Bretschneider spectrum, which belongs to a two parameter spectrum family that permits 
period and wave height to be assigned separately, should be considered for future 


research. 
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APPENDIX: EIGHTH-ORDER FILTER 


function [rho] = predict (delta) 

Matlab script to function as a Eighth-Order Digital Filter for 
Predicting Future Seaway Elevation Response 
Tracks and matches Pierson Moskowitz/Pressure Profile 
Spectrum and predicts responses one full wave length 


Variables 

Fignificant wave height 

Frequency 

Defines Pierson Moskowitz expression for fully developed seas 
Time vector 
a = Damping ratio 

Period interval 
A,B = Continuous plant model 

Phi,Gamma = Discrete plant model 

K = Filter Gains 

error= 

x2 = Estimate of state vector 

yY = System’s output 

rho = Cross-correlation coefficient 

delta is the number of time steps, (i.e. Phi is delta minus one) 


an rMNes 
rou wile 


de dP dP dP dP dP dO dP dO de dP dP dP dP de dP dP dP dP dP 


h=3; 
w=[0.3:0.05:3];dw=0.05; 
fl,msj=size(w); 
S=w; 
for i=i:ms 
S(i)=8.1e-3*32.2°2/ (w(i) *5) *exp(-33.56/h*2/ (w(i)%4)); 
end; 
ws={0.3:0.05:3); 


lambda=(171.61,146.52,127.62,112.86,100.98,91.21,83.01,76.02,70,64.73,... 
60.08,55.94,52.23,48.88,45.83,43.05,40.49,38.14,35.96,33.93,32.05,30.29,... 
28.64.27.1,25.66,24.3,23.03,21.84,20.71,19.66,18.67,17.75,16.88,16.07,15.30,... 
14.$9,13.92,13.29,12.7,12.15,11.63,11.15,10.69,10.26,9.86,9.47,9.41,8.77,... 
8.45,8.15,7.86,7.59,7.33,7.08,6.85} ; 


% Pressue Measurement Simulation 


Llambda=lambda.*3.28; 
t=(0:0.1:200]; 
Y=zeros(l1,length(t)); 
for i=l:length(ws) ; 
phi=rand(1,length(ws) ) ;phi=phi-mean (phi) ; 
y(i,:)=(cosh(6*pi/lambda (i) )/cosh(40*pi/lambda(i)))... 
* (sqrt (S(i) *2*dw) )*cos (ws (i) *t+phi(i)*pi*2); 
end; 
for j=l:length(ws) ; 
Y=H=Y¥+y(J,:); 
end; 


% Defines Coefficients of transfer function of Eighth-Order Filter 


T=0.1;w0=1.28;zeta=0.4; 


n=[1/w0,0];d=[(1/w0)*%2,2*zeta/w0,1}; % Defines Fourth-Order Expression 
num=conv (n,n) ;num4=conv (num,num) ; 
den=conv (d,d) ;den4=conv (den, den) ; % Defines transfer function 
% Defines innovater gains for subject filter 
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{A1l,B1,C1,Dlj=tf2ss (num4,den4) ; 

([Phi,Gamma]=c2d(A1,B1,T); 

K=@lqr(Phi’,C1’, eye(8,8)*10,.1); eigl=abs(eig(Phi-K’*C1)) 
x2=zeros(8,length(t));x4=x2;P2=zeros(i,length(t));P6=P2; 


% For loop for closed loop filter 


for 1=1L>Llengtate) 
x2(:,itl)=(Phi) *x2(:,1)+K’* (¥(1)-C1*x2(:,1)); 
x60(:,i+1)=(Phi*delta) *x2(:,1+1); 
P2(i)=C1*x2(:,1); 
P60(i)=C1*x60(:,1); 

end; 


for i=1:{length(t)-delta); 
error60(i)=(P60(i)-Y(itdelta));Y60(i)=Y(i+delta) ; 


end; 
% Calculates cross-correlation coefficient between P-M spectrum and filte 


bo 


rho= (Y60*P60(1: (2001-delta) )‘’/(std(Y60)*std(P60)))/(2001-delta) ; 
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